// -----------------------------------------------------------------------------
// M A R I T I M E   R E S E A R C H   I N S T I T U T E   N E T H E R L A N D S
// -----------------------------------------------------------------------------
// Program  : lin_wacu
// File     : lin_wacu.cpp
// -----------------------------------------------------------------------------
// Revisions:
// Februari 13, 2009  asicherer
// -----------------------------------------------------------------------------
// testbed for the interactive waves model
// -----------------------------------------------------------------------------

// note: as this is testbed-only code: the architecture is not to be taken as an example

#define REMAIN // leave display on after simulation

#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <deque>
#include <vector>
//#include "FloatUtils.h"
#include "Array2D.h"
#include<math.h>
#include "defines.h"
#include "SolverFactory.h"
#include "RrbAkashSolver.h"
#include "RrbSolver.h"
#include <sstream>
#define PI 3.14159265
#define Chorin
#define openmp2

static REAL currenttime = 0.0; ///< time stamp
static REAL timestep; ///< duration of one step in seconds
static int Nx1, Nx2; ///< size in points in x1,x2 direction
static REAL dx1, dx2; ///< grid step size in x1,x2 direction
static REAL grid; ///< minimum of grid step sizes
static REAL g; ///< gravity
static bool running = true; // simulation running?

static RrbAkashSolver *RrbSolver_Akash=NULL;

/// reads test file
static void setup(Array2D<REAL> *pHeight , Array2D<REAL> *pSpeed , const REAL &dt, const int &flag)
{

    /* -------------------------------------------------------------------*/
    // The System Data
  /*  Nx1 = 400;
    Nx2 = 250;
    REAL  Lx1 = 200;
    REAL Lx2 = 125;*/


    Nx1 = 20;
    Nx2 = 20;
    REAL  Lx1 = 100;
    REAL Lx2 = 100;
//

    REAL Ux = 1.0 ;
    REAL Uy =1.0;
    REAL beta = 0.5;
    REAL h=50;
    REAL time_max =80;

    g =9.81;
    REAL tolerance = 0.000001;
    timestep=dt;
    int steps =int(time_max/timestep);
    //steps=10;
    int COLUMNS = Nx1+3;
    int ROWS = Nx2+3;
    dx1 = Lx1 / (Nx1-1);
    dx2 = Lx2 / (Nx2-1);
    REAL dt1=1.0/timestep;
    std::cout << dt1 << std::endl;
    REAL dt_cfl = std::min(dx1,dx2)*0.5/sqrt(g*h);
    std::cout << "Delta t CFL is:" << dt_cfl << std::endl;
    std::cout << "Delta t Ratio is:" << timestep/dt_cfl << std::endl;



    // Construct the System Matrices

    Array2D<REAL> mat_zz_c(COLUMNS, ROWS) ;
    Array2D<REAL> mat_zz_w(COLUMNS, ROWS);
    Array2D<REAL> mat_zz_s(COLUMNS, ROWS);

    Array2D<REAL> mat_ff_c(COLUMNS, ROWS);
    Array2D<REAL> mat_ff_w(COLUMNS, ROWS);
    Array2D<REAL> mat_ff_s(COLUMNS, ROWS);

    Array2D<REAL> mat_zf_c(COLUMNS, ROWS);
    Array2D<REAL> mat_zf_w(COLUMNS, ROWS);
    Array2D<REAL> mat_zf_s(COLUMNS, ROWS);

    Array2D<REAL> mat_fz_c(COLUMNS, ROWS);
    Array2D<REAL> mat_fz_w(COLUMNS, ROWS);
    Array2D<REAL> mat_fz_s(COLUMNS, ROWS);

    mat_zz_c =1.0;
    mat_zz_w =0.0;
    mat_zz_s =0.0;

    mat_ff_c =1.0;
    mat_ff_w =0.0;
    mat_ff_s =0.0;

    mat_zf_c =1.0;
    mat_zf_w =0.0;
    mat_zf_s =0.0;

    mat_fz_c =1.0;
    mat_fz_w =0.0;
    mat_fz_s =0.0;



    for (int i2 = 1; i2 <= Nx2; ++i2)
    {
        for (int i1 = 1; i1 <= Nx1; ++i1)
        {
            mat_zz_c[i1][i2] = dt1;
            mat_zz_w[i1][i2] = -beta*Ux*0.5/dx1;
            mat_zz_s[i1][i2] = -beta*Uy*0.5/dx2;


            mat_ff_c[i1][i2] = dt1;
            mat_ff_w[i1][i2] = -beta*Ux*0.5/dx1;
            mat_ff_s[i1][i2] = -beta*Uy*0.5/dx2;

            mat_zf_c[i1][i2] = -beta*h*(2/(dx1*dx1)+2/(dx2*dx2));
            mat_zf_w[i1][i2] = beta*h/(dx1*dx1);
            mat_zf_s[i1][i2] = beta*h/(dx2*dx2);

            mat_fz_c[i1][i2] = beta*g;
            mat_fz_w[i1][i2] = 0;
            mat_fz_s[i1][i2] = 0;

        }
    }

    // Apply the Boundary Condition- Neumann
    /* north and south */
    for (int i1 = 1; i1 <= Nx1; ++i1)
    {
        mat_zz_c[i1][Nx2] -=mat_zz_s[i1][Nx2];
        //mat_zz_cn[i1][Nx2]=0;
        mat_zz_c[i1][1] +=mat_zz_s[i1][1];
        mat_zz_s[i1][1]=0;

        mat_ff_c[i1][Nx2] -=mat_ff_s[i1][Nx2];
        //  mat_ff_n[i1][Nx2]=0;
        mat_ff_c[i1][1] +=mat_ff_s[i1][1];
        mat_ff_s[i1][1]=0;

        mat_zf_c[i1][Nx2] +=mat_zf_s[i1][Nx2];
        // mat_zf_n[i1][Nx2]=0;
        mat_zf_c[i1][1] +=mat_zf_s[i1][1];
        mat_zf_s[i1][1]=0;

        mat_fz_c[i1][Nx2] +=mat_fz_s[i1][Nx2];
        // mat_fz_n[i1][Nx2]=0;
        mat_fz_c[i1][1] +=mat_fz_s[i1][1];
        mat_fz_s[i1][1]=0;
    }

    /* east and west*/
    for (int i2 = 1; i2 <= Nx2; ++i2)
    {
        mat_zz_c[Nx1][i2] -=mat_zz_w[Nx1][i2];
        //  mat_zz_e[Nx1][i2]=0;
        mat_zz_c[1][i2] +=mat_zz_w[1][i2];
        mat_zz_w[1][i2]=0;

        mat_ff_c[Nx1][i2] -=mat_ff_w[Nx1][i2];
        //  mat_ff_e[Nx1][i2]=0;
        mat_ff_c[1][i2] +=mat_ff_w[1][i2];
        mat_ff_w[1][i2]=0;

        mat_zf_c[Nx1][i2] +=mat_zf_w[Nx1][i2];
        //  mat_zf_e[Nx1][i2]=0;
        mat_zf_c[1][i2] +=mat_zf_w[1][i2];
        mat_zf_w[1][i2]=0;

        mat_fz_c[Nx1][i2] +=mat_fz_w[Nx1][i2];
        //  mat_fz_e[Nx1][i2]=0;
        mat_fz_c[1][i2] +=mat_fz_w[1][i2];
        mat_fz_w[1][i2]=0;
    }


    RrbSolver_Akash = new RrbAkashSolver(mat_zz_c , mat_zz_s, mat_zz_w,
                                   mat_ff_c , mat_ff_s, mat_ff_w,
                                   mat_zf_c , mat_zf_s, mat_zf_w,
                                   mat_fz_c , mat_fz_s, mat_fz_w,
                                   Nx1, Nx2, tolerance);

    // Initiate the Variable Zeta and Phi..
    // Allocate memory to pHeight and pSpeed

    pHeight =   new Array2D<REAL> (COLUMNS, ROWS);
    pSpeed =    new Array2D<REAL> (COLUMNS, ROWS);

    Array2D<REAL> &m_Zeta = *pHeight;
    Array2D<REAL> &m_Phi = *pSpeed;
    Array2D<REAL> m_Psi(COLUMNS, ROWS);
    m_Zeta =0;
    m_Phi=1;
    m_Psi=1;
    for (int i1 = 0; i1 <= Nx1+2; ++i1)
    {
        for (int i2 = 0; i2 <= Nx2+2; ++i2)
        {
            m_Zeta[i1][i2] = cos(4*PI*(i1-1)*dx1/Lx1);
            m_Phi[i1][i2] = 1;
            m_Psi[i1][i2] = 1;

        }
    }
    //-------------------------------------------------------------------------------------------------

    // Construct the RHS Stencil
    REAL beta1 =1.0-beta;

    for (int i2 = 1; i2 <= Nx2; ++i2)
    {
        for (int i1 = 1; i1 <= Nx1; ++i1)
        {
            mat_zz_c[i1][i2] = dt1;
            mat_zz_w[i1][i2] = beta1*Ux*0.5/dx1;
            mat_zz_s[i1][i2] = beta1*Uy*0.5/dx2;


            mat_ff_c[i1][i2] = dt1;
            mat_ff_w[i1][i2] = beta1*Ux*0.5/dx1;
            mat_ff_s[i1][i2] = beta1*Uy*0.5/dx2;

            mat_zf_c[i1][i2] = beta1*h*(2/(dx1*dx1)+2/(dx2*dx2));
            mat_zf_w[i1][i2] = -beta1*h/(dx1*dx1);
            mat_zf_s[i1][i2] = -beta1*h/(dx2*dx2);

            mat_fz_c[i1][i2] = -beta1*g;
            mat_fz_w[i1][i2] = 0;
            mat_fz_s[i1][i2] = 0;
        }
    }
    // Apply The boundary Condition on the Stencil!
    for (int i1 = 1; i1 <= Nx1; ++i1)
    {
        mat_zz_c[i1][Nx2] -=mat_zz_s[i1][Nx2];
        //mat_zz_cn[i1][Nx2]=0;
        mat_zz_c[i1][1] +=mat_zz_s[i1][1];
        mat_zz_s[i1][1]=0;

        mat_ff_c[i1][Nx2] -=mat_ff_s[i1][Nx2];
        //  mat_ff_n[i1][Nx2]=0;
        mat_ff_c[i1][1] +=mat_ff_s[i1][1];
        mat_ff_s[i1][1]=0;

        mat_zf_c[i1][Nx2] +=mat_zf_s[i1][Nx2];
        // mat_zf_n[i1][Nx2]=0;
        mat_zf_c[i1][1] +=mat_zf_s[i1][1];
        mat_zf_s[i1][1]=0;

        mat_fz_c[i1][Nx2] +=mat_fz_s[i1][Nx2];
        // mat_fz_n[i1][Nx2]=0;
        mat_fz_c[i1][1] +=mat_fz_s[i1][1];
        mat_fz_s[i1][1]=0;
    }

    /* east and west*/
    for (int i2 = 1; i2 <= Nx2; ++i2)
    {
        mat_zz_c[Nx1][i2] -=mat_zz_w[Nx1][i2];
        //  mat_zz_e[Nx1][i2]=0;
        mat_zz_c[1][i2] +=mat_zz_w[1][i2];
        mat_zz_w[1][i2]=0;

        mat_ff_c[Nx1][i2] -=mat_ff_w[Nx1][i2];
        //  mat_ff_e[Nx1][i2]=0;
        mat_ff_c[1][i2] +=mat_ff_w[1][i2];
        mat_ff_w[1][i2]=0;

        mat_zf_c[Nx1][i2] +=mat_zf_w[Nx1][i2];
        //  mat_zf_e[Nx1][i2]=0;
        mat_zf_c[1][i2] +=mat_zf_w[1][i2];
        mat_zf_w[1][i2]=0;

        mat_fz_c[Nx1][i2] +=mat_fz_w[Nx1][i2];
        //  mat_fz_e[Nx1][i2]=0;
        mat_fz_c[1][i2] +=mat_fz_w[1][i2];
        mat_fz_w[1][i2]=0;
    }


    Array2D<REAL> cn = mat_zz_s;
    Array2D<REAL> ce = mat_zz_w;

    for (int jj = 0; jj <= Nx2 + 1; ++jj)
    {
        for (int ii = 0; ii <= Nx1 + 1; ++ii)
        {
            cn[ii][jj] = -mat_zz_s[ii ][jj + 1]; // Make use of the fact that the matrix is Skew Symmetric.
            ce[ii][jj] = -mat_zz_w[ii + 1][jj ]; // Make use of the fact that the matrix is Skew Symmetric.
        }
    }

    std::vector<Array2D<REAL> > mat_zz;
    mat_zz.push_back(mat_zz_c);
    mat_zz.push_back(mat_zz_s);
    mat_zz.push_back(mat_zz_w);
    mat_zz.push_back(cn);
    mat_zz.push_back(ce);

    cn = mat_ff_s;
    ce = mat_ff_w;

    for (int jj = 0; jj <= Nx2 + 1; ++jj)
    {
        for (int ii = 0; ii <= Nx1 + 1; ++ii)
        {
            cn[ii][jj] = -mat_ff_s[ii ][jj + 1]; // Make use of the fact that the matrix is Symmetric.
            ce[ii][jj] = -mat_ff_w[ii + 1][jj ]; // Make use of the fact that the matrix is Symmetric.
        }
    }

    std::vector<Array2D<REAL> > mat_ff;
    mat_ff.push_back(mat_ff_c);
    mat_ff.push_back(mat_ff_s);
    mat_ff.push_back(mat_ff_w);
    mat_ff.push_back(cn);
    mat_ff.push_back(ce);

    cn = mat_zf_s;
    ce = mat_zf_w;

    for (int jj = 0; jj <= Nx2 + 1; ++jj)
    {
        for (int ii = 0; ii <= Nx1 + 1; ++ii)
        {
            cn[ii][jj] = mat_zf_s[ii ][jj + 1]; // Make use of the fact that the matrix is Symmetric.
            ce[ii][jj] = mat_zf_w[ii + 1][jj ]; // Make use of the fact that the matrix is Symmetric.
        }
    }
    std::vector<Array2D<REAL> > mat_zf;
    mat_zf.push_back(mat_zf_c);
    mat_zf.push_back(mat_zf_s);
    mat_zf.push_back(mat_zf_w);
    mat_zf.push_back(cn);
    mat_zf.push_back(ce);


    cn = mat_fz_s;
    ce = mat_fz_w;

    for (int jj = 0; jj <= Nx2 + 1; ++jj)
    {
        for (int ii = 0; ii <= Nx1 + 1; ++ii)
        {
            cn[ii][jj] = mat_fz_s[ii ][jj + 1]; // Make use of the fact that the matrix is Symmetric.
            ce[ii][jj] = mat_fz_w[ii + 1][jj ]; // Make use of the fact that the matrix is Symmetric.
        }
    }

    std::vector<Array2D<REAL> > mat_fz;
    mat_fz.push_back(mat_fz_c);
    mat_fz.push_back(mat_fz_s);
    mat_fz.push_back(mat_fz_w);
    mat_fz.push_back(cn);
    mat_fz.push_back(ce);

    std::cout << "entered the time loop" << std::endl;


    // Construct the Stencil for the matrix S_ZetaPsi
    std::vector<Array2D<REAL> > mat_zp;
    REAL coeff_x = h*h/(3*dx1*dx1);
    REAL coeff_y = h*h/(3*dx2*dx2);
    REAL coeff_center = -2*(coeff_x + coeff_y);
    constructMatrix(Nx1, Nx2, coeff_x, coeff_y, coeff_center ,mat_zp );


    // Construct the Stencil for the matrix S_PsiPsi
    std::vector<Array2D<REAL> > mat_pp;
    coeff_x = -2*h*h*h/(15*dx1*dx1);
    coeff_y = -2*h*h*h/(15*dx2*dx2);
    coeff_center = -2*(coeff_x + coeff_y) + h/3;
    constructMatrix(Nx1, Nx2, coeff_x, coeff_y, coeff_center ,mat_pp );

    Array2D<REAL> mat_pp_c = mat_pp[0];
    Array2D<REAL> mat_pp_s = mat_pp[1];
    Array2D<REAL> mat_pp_w = mat_pp[2];


    // The matrix for Stencil S_PsiPhi is same as S_ZetaPsi
    std::vector<Array2D<REAL> > mat_pf(mat_zp);
    //RrbPsi = new RrbSolver(&mat_pp_c,&mat_pp_s, &mat_pp_w, Nx1, Nx2, tolerance);

    Solver *RrbPsiSolver = SolverFactory::create("rrb", &mat_pp_c,&mat_pp_s, &mat_pp_w, Nx1, Nx2, tolerance);


    // -----------------------------------------------------------------------------------

    // Start the Time Loop
    for (int t1 =1 ; t1<= steps ; ++ t1)
    {
        // Apply Boundary Conditions
      //  neumannBoundary(&m_Zeta,Nx1, Nx2);
      //  neumannBoundary(&m_Phi,Nx1,Nx2);


        // Compute RHS
        Array2D<REAL> RHS_ZZ = MatAdd(matvec(mat_zz,m_Zeta), matvec(mat_zf,m_Phi));
        Array2D<REAL> RHS_FF = MatAdd(matvec(mat_fz,m_Zeta), matvec(mat_ff,m_Phi));

        double wavetime = seconds();

        // The Solver
        if (flag==1)
        {

            RrbSolver_Akash->solve(&RHS_ZZ,&RHS_FF,&m_Zeta,&m_Phi);
        }
        else
        {
            RrbSolver_Akash->GMRES(&RHS_ZZ,&RHS_FF,&m_Zeta,&m_Phi);
        }
        double wavetime1= seconds();
        // Reapply the Boundary condition on Phi
     //   neumannBoundary(&m_Phi,Nx1,Nx2);
#ifdef Chorin
        Array2D<REAL> RHS_PP = scalarmult(-1,matvec(mat_pf, m_Phi));
        // Apply the Boundary Condition on Psi
      //  neumannBoundary(&m_Psi,Nx1,Nx2);

        // Solve the Elliptic System
        RrbPsiSolver->solve(&RHS_PP,&m_Psi);
        double wavetime2 = seconds();

        // Correct Zeta
        m_Zeta = MatAdd(m_Zeta,scalarmult(timestep,matvec(mat_zp,m_Psi)));
#endif

      //  std::cout << "Time for Solver" << flag << " : " << wavetime1-wavetime << "Time for CG" << wavetime2-wavetime1 <<std::endl;
      //  std::ostringstream ss;
     //   ss << t1;
      //  std::string s = "z" + ss.str();
      // // PrintArray2Dcontents1D(s,m_Zeta,1,Nx1,1,Nx2);
    }
    PrintArray2Dcontents("",m_Zeta,1,Nx1,1,Nx2);
}




/// test program
int main(int argc, char *argv[])
{

    // Need to define pHeight as a pointer, because we do not yet know the problem size.
    Array2D<REAL> *pHeight=0;
    Array2D<REAL> *pSpeed=0;


    // Check Matrix operations:

    /*Array2D<REAL> A1 (10,10);
    Array2D<REAL> A2 (10,10);

    A1=3;
    A2=2;

    REAL s = dotproduct(A1, A2);
    std::cout << s << std::endl;
    A2[1][1]=1.0;
    Array2D<REAL> A3 = MatSub(A1,A2);
    A3 = scalarmult(5,A3);
    PrintArray2Dcontents1D("",A3,0,A3.dim1()-1,0,A3.dim2()-1);*/


    REAL dt=0.01;
    if (argc >=1)
         dt = atof(argv[1]);
    std::cout << dt << std::endl;
    int flag=1;
    //if (argc >= 2)
    flag = atoi(argv[2]);

    double waveTime =seconds();
    setup(pHeight, pSpeed,dt, flag);
    waveTime = seconds() - waveTime;
    std::cout << "Total time is" << waveTime << std::endl;

    return 0;

}

